Nonequilibrium Response from the dissipative 
Liouville Equation 

Matteo Colangeli| 

School of Mathematical Sciences, Queen Mary University of London, Mile End Road, 
London El 4NS, UK 

E-mail: colangeliOcalvino . polito . it 

Abstract. The problem of response of nonequilibrium systems is currently under 
intense investigation. We propose a general method of solution of the Liouville 
Equation for thermostatted particle systems subjected to external forces which retains 
only the slow degrees of freedom, by projecting out the majority of fast variables. 
Response formulae, extending the Green-Kubo relations to dissipative dynamics are 
provided, and comparison with numerical data is presented. 

PACS numbers: 05.20.Dd,05.20.Gg, 05.70.Ln 



| Present address: Dipartimento di Matematica, Politecnico di Torino, Corso Duca degli Abruzzi 24, 
1-10129 Torino, Italy 



Nonequilibrium Response from the dissipative Liouville Equation 



2 



1. Introduction 



The celebrated theory of linear response of systems close to equilibrium has been 
developed in the first half of the XX century, and has proven highly successful in 
describing a very wide range of macroscopic phenomena, particularly thanks to the 
usefulness of the Fluctuation-Dissipation Theorems (FDT) [1]. In the second half of the 
XX century, numerous generalizations have been proposed to extend the classical theory 
beyond the linear regime, by considering the microscopic evolution as represented by 
deterministic dynamical systems, by dynamical systems perturbed by some appropriate 
noise, or by purely stochastic processes, which reach a nonequilibrium steady state. 
In particular, relations known as Fluctuation Relations have been introduced together 
with different nonequilibrium generalizations of the FDT. One then investigates the 
variations of these dynamics under certain perturbations. Remarkably, it turns out that 
rather different approaches lead to similar results, cf . Ref . [2] for recent reviews of the 
historical developments and of the recent results. The recent developments, however, do 
not constitute a consistent and comprehensive theory yet, therefore a variety of different 
approaches is worthy investigating. For instance, Ruelle has shown [3] that, for uniformly 
hyperbolic dynamical systems, the linear response away from equilibrium is very similar 
to the linear response close to equilibrium: the Kramers-Kronig dispersion relations hold, 
and the FDT survives in a modified form, which accounts for the oscillations around 
the relevant attractors. If the chaotic hypothesis does not hold, Ruelle concludes that 
two new phenomena may arise. The first is a violation of linear response, in the sense 
that the nonequilibrium steady state does not depend differentiably on parameters. 
The second phenomenon is a violation of the dispersion relations: the susceptibility 
has singularities in the upper half complex plane. These 'acausal' singularities are 
actually due to 'energy nonconservation': for a small periodic perturbation of the 
system, the amplitude of the linear response is arbitrarily large. This means that the 
steady state of the dynamical system under study is not 'inert' but can give energy 
to the outside world, something rather different from the behavior of an equilibrium 
state. This approach is based on the smoothness property of SRB measures [4] along 
the unstable fibres of the attractors, something hardly realizable in the phase space 
probability distribution of physically relevant nonequilibrium systems. Therefore, this 
approach is limited to the properties of systems whose departure from the ideal uniform 
hyperbolicity does not seriously affect the observed behavior [3]. From a different 
perspective, Vulpiani et al. [5] have studied the response of dynamical systems to finite 
amplitude perturbation, assuming that their phase space probability distributions are 
not singular with respect to the Lebesgue measure, hence have an integrable density. 
These authors found that a generalized fluctuation-response relation holds, and it links 
the average relaxation towards the steady state to the invariant measure of the system. 
Furthermore, this approach points out the relevance of the amplitude of the initial 
perturbation, something which does not pertain to the approach mentioned above. The 
use of regular distributions is justified by considering that any physical phenomenon 
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is affected by some noise, capable of smoothing out the singularities produced by the 
dissipative deterministic dynamics. However, one may also like to consider the noise 
as emerging from a projection procedure [6], or may investigate the possibility that the 
distribution (not necessarily invariant) be expressed in terms of appropriate expansions 
around a reference one. In this work we pursue the latter approach and propose a 
method which resembles the Chapman- Enskog technique [7] for the Boltzmann Equation 
(BE), which is suitable for the generalized Liouville Equation (LE) of deterministically 
thermostatted particle systems. The method takes inspiration from the Bogolyubov 
hypothesis [8] of time scales separation and postulates that the dynamics of the single- 
particle distribution function is driven by the dynamics of a set of selected fields [9]. 
Analogously, in the context of many-particle systems, one can conjecture that the 
dynamics of the phase space probability density is triggered by the dynamics of some 
averaged phase space quantities. In particular, in the theory of the BE, the notion of 
slow variables was corroborated by the assumption of local equilibrium. At the level of 
the LE, on the other hand, a similar notion can be introduced under the assumption of 
a weakly dissipative dynamics, in order to let the resulting probability density be not 
too far from the Gibbsian density, which is required for the application of the method. 

2. The Chapman-Enskog method of solving the Boltzmann Equation 

In this Section we review some of the main features of the Chapman-Enskog method, 
which we intend, next, to employ, in Sec. 3, in the context of thermostatted particle 
systems, to study the evolution of the probability density in the full phase space. The 
approach, introduced by Enskog and popularized in a modified version by Chapman and 
Cowling [7], consists of a mathematical procedure meant to approximate the solution 
of the BE. The method deserves a special mention in kinetic theory of gases since the 
Boltzmann pioneering works, as it allows to derive the Navier-Stokes-Fourier equations 
of hydrodynamics from the BE, by elucidating and making use of concepts, such as 
the scaling (hydrodynamic) limit of a kinetic equation, the notion of local equilibrium 
and the time scales separation, which led to the systematic construction of projection 
operator methods in statistical mechanics as well as to the later development of special 
large deviation methods in the theory of stochastic processes [10]. Let us start by writing 
down the BE, which (in absence of external forces) reads: 

dtf = -v • V/ + C(f) (1) 

where / = /(r, v, t) is the one particle distribution function, depending on position 
r, particle velocity v and time t and C(f) is a nonlinear integral collision operator. By 
rescaling the BE with characteristic time and length scales, the Knudsen number (- = \ 
arises naturally from (1), which, in dimensionless form, becomes: 

d t f = -v ■ V/ + -C(f) (2) 
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where L denotes a macroscopic length scale (e.g. the sice of the system) and A is the 
mean free path, i.e. the average distance covered by a moving particle between successive 
collisions with other moving particles. In the limit e — >■ 0, which is commonly referred to 
as the hydrodynamic limit, the fluid becomes dense enough that the dynamics, described 
by (2), is dominated by the effect of the collisions. By introducing the notion of local 
equilibrium, (i.e. the existence of macroscopic regions each of which is in an equilibrium 
state which may be different from the equilibrium states found in other regions) it 
is possible to define, locally, a set of selected variables x, known as hydrodynamic 
fields (which correspond to the collision invariants: number of particles density n(r, t), 
momentum density n(r, t)u(r, t), and kinetic energy density e(r,t) = |ro(r, t)/sgT(r, t)), 
which enter the definition of Maxwell-Boltzmann statistics: 

fn o m(V-U(r,t)) 2 

where m is the mass of the particle. The local Maxwellians (3) are not, in general, 
solutions of (1). Nevertheless, in the hydrodynamic limit, where the length scale of the 
spacial inhomogeneities, L ~ e _1 , tends to diverge, they offer a good approximation to 
the exact solution. In fact, it is evident from Eq. (2) that for e — > the last term on the 
right becomes singular. Then, the only way to avoid the singularity is that the collision 
term itself vanishes, which is guaranteed by the form of the local Maxwellians. It is 
worth to notice that Eq. (2) can also be considered, in the Fourier-Laplace space, as an 
eigenvalue problem for the operator A = C(f) — zek • v. As discussed in Ref. [11], some 
technical problems arise when one attempts to solve the eigenvalue problem away from 
the strict hydrodynamic limit e — > 0, and alternative techniques, based on projection 
operators, have been recently made available [12]. However, in those circumstances 
where e <C 1, it makes sense to attempt a perturbative method to solve Eq. (2), by 
introducing the following expansion: 

oo 

/ = £e7 W ( 4 ) 

1=0 

Thus, by inserting (4) into (2), one may define a microscopic time derivative Q™ %cro f as: 



d't 



oo 1 oo 

<*°f = -v • v[J2 f (l) A + -c(J2 / ( ° e ') ( 5 ) 



1=0 1=0 



Following [12], in the theory of the BE it is also possible to introduce a macroscopic 
time derivative, d" iacro f, provided that the assumption of local equilibrium holds. 
The idea is to employ the notion of normal solutions, which allows to express the 
spatial and temporal dependence of the terms f^ through the hydrodynamic fields, 
f^(r,v,t) = /^(x(r, t), v). By applying a sort of chain rule, we can, then, write: 

d(Y°°_ e j f (i) ) 

grnacroj = Q t X (6) 
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The time derivative of the hydrodynamics fields, in (6), is found to be given by: 

oo 

a,x = ^ e fe 9f ) x (?) 

fc=0 

where d\ ! x = — J w(v)v • V/^oPv and where the w's are lower order Sonine 
polynomials of v. The method, then, requires to equalize the two time derivatives 
Qmicro j an( j Qmacro j (supplemented by (7)). The equality, in fact, provides a cascade of 
equations for the various terms of the expansion /^(x(r, t), v), whose leading order is 
represented by the family of Maxwellians (3), and whose first correction, can be 
obtained [7] from: 



j w(v)v • Wf {Q \v)d\ = v • V/ (0) - (8) 



The lowest order of the equations of hydrodynamics (7) is known as Euler 
hydrodynamics and features vanishing transport coefficients, whereas, by adding 
the first correction one obtains the dissipative Navier-Stokes- Fourier equations, 
which are endowed with constitutive expressions for the stress tensor and the 
heat flux [12]. However, higher order corrections of the Chapman- Enskog method, 
resulting in hydrodynamic equations with higher derivatives (Burnett and super-Burnett 
hydrodynamic equations) face severe difficulties both from theoretical, as well as from 
the practical point of view and various regularization methods have been suggested 
[9, 13]. 

3. Dissipative Liouville Equation with external forcing and nonequilibrium 
response 

The method discussed in the previous Section and applied in the context of kinetic 
theory can be extended to Hamiltonian particle systems, as it is formally possible to 
construct perturbative techniques to solve the LE. Such perturbation theories allow 
to compute the higher order corrections to the phase space counterpart of the local 
equilibrium in the BE and lead to expressions of transport coefficients in terms of the 
detailed microscopic dynamics. The standard method which succeeds in this derivation 
was originally proposed by Green and Kubo [14]. Any perturbation theory based on an 
expansion of the probability density in powers of a small parameter e needs to extend and 
justify the assumption of local equilibrium in phase space and to prove that higher order 
corrections are suitable to refine the lower order approximations. This is a delicate 
issue, as, formally, there is no indication of the convergence of such an expansion 
in phase space and, mostly, there is no analog of the Caflish Theorem [15], formally 
proving that a proper truncation of the series (4) approximates the solution to the BE. 
Furthermore, as pointed out by the Authors in [16], in order to obtain finite transport 
coefficients, one needs to postulate a strong decay property of the time autocorrelations 
of mass, momentum and energy currents. We intend to show, here, that, in the context 
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of externally driven thermostatted particles systems, the analogy with kinetic theory 
can be, to some extent, retained and most of the concepts formerly employed in the 
derivation of hydrodynamics from the BE, can be used to obtain nonequilibrium response 
formulae from the LE. In particular, the dimensionless parameter e which arises after a 
proper rescaling of the LE, does not constitute just a mathematical book-keeping device 
which is eventually set to unity, but it can be endowed with a physical content, such as 
it was in the case of the BE. We consider a particle system, held, for time t G (— oo, 0], 
in equilibrium with an external bath at temperature T. At time t = an external force 
F ea:t , starts acting and, hence, some energy is "pumped" inside the system. At the same 
time a deterministic thermostat is switched on in order to remove part (or the whole) 
of the energy provided by the external force, and to achieve, on a large time scale, 
a nonequilibrium steady state. The dynamics we have in mind are the deterministic 
thermostatted dynamics of nonequilibrium molecular dynamics [17]. In particular, we 
refer to those systems for which the dimensionless LE takes the form: 

p = -f ■ V/9 + enp (9) 

with: 

i = f o + eR(F) = S • V# + eR(T) (10) 

where S is the symplectic matrix, Ho is the Hamiltonian of the many-particle system, 
r denotes the corresponding Hamiltonian contribution to the dynamics and R(T) is the 
term which spoils the conservativity of the equations of motion. Furthermore, in (9), 
we introduced a shorthand notation for the phase space contraction rate: k = — V ■ F 
and e is a dimensionless parameter which is proportional to the intensity of the external 
driving. A paradigmatic example, where the parameter e attains a sensible structure, is 
provided by the Gaussian thermostatted Lorentz gas of hard spheres, with Hamiltonian 
Hq(F) = y~. By adding an external force, here represented by F ext = gE, with q the 
electric charge and E the external electric field, and the thermostat, both entering the 
definition of R, in (10), one finds: R = gE — a(p)p, where m the mass of the particle and 
a = q p2 ^ = k{F). For simplicity, we consider just a single particle and we also do not 
take the presence of scatterers explicitly into account, as it does not contribute to the 
phase space contraction rate, which, as discussed below, is the only observable we are 
going to be concerned with, in our model. By rescaling all dimensional quantities with 
proper characteristic scales, the resulting dimensionless LE takes precisely the structure 
indicated in (9), and, in analogy with the definition of the Knudsen number introduced 
in Sec. 2, the parameter e is, then, given by the ratio of the energy absorbed by the 
thermostat, A£d iss , to a characteristic energy of the system K: 

e= 11 
it K 

where denotes E a typical intensity of the applied electric field, L a reference 
length scale, v a characteristic velocity of the dynamics and K may be identified by the 
kinetic energy of the system or, in systems enjoying local thermodynamic equilibrium, 
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may be also related to the temperature T of the system. We observe that in the 
Zubarev's seminal papers [18], an infinitesimal source term, corresponding to our e, 
was introduced ad hoc in the construction of the nonequilibrium statistical operator 
from the LE, whereas, here, in the context of thermostatted particle systems, a small 
parameter arises naturally after the adimensionalization. The parameter e can be small 
or large, depending on the details of the coupling with the thermostat. It is well known 
that, in strongly dissipative systems (e ^> 1), the support of the invariant measure is an 
attractor endowed with a very thin fractal structure, which appears strongly at variance 
with a regular distribution. On the other hand, in the limit of weak dissipation (i.e.: 
weak coupling with the thermostat), e <C 1, it makes sense to attempt a perturbative 
technique to solve Eq. (9), qhich is performed by expanding p in powers of e: 

oo 

p=£ e y° (12) 

1=0 

where p^ obeys a purely conservative equilibrium dynamics (unaffected by the external 
field and by the thermostat): 

p(°)(t) = e - A V 0) (0) (13) 

with A = To • V the Liouvillian operator. Let us point out, here, the analogy with Sec. 
2. In the context of the BE, we have already remarked that, in the limit e — > oo, it holds 
/ — > f , i.e., the local equilibria (3) are a good approximation of the exact statistics 
of the system. Similarly, we find, here, that for vanishing e, p — > p(°\ i.e. the density 
reduces to the purely Hamiltonian contribution. 
With (10) and (12), Eq. (9) transforms into: 

p = _r . V(p (0) + ep« + ...) - e{R ■ V - ^)(p (0) + ep« + ...) (14) 

Moreover, if the assumption of time scales separation holds, one may wish to 
decompose the dynamics of the probability density into a contribution (expressed via a 
projection operator P, to be defined below in (16)) which is solely dependent on some 
relevant slow variables, plus a contribution given in terms of an orthogonal projector 
Q = 1 — P (with 1 the identity operator), which features the dynamics of the fast 
variables. The decomposition of the dynamics of the density as induced by a specific 
choice of the projector P was provided by Robertson [20]: 

p = Pp + [ K(t,t')p(t')dt' (15) 
Jo 

where the term K(t,t')p(t') includes the contribution of fluctuations of the 
overwhelming majority of the fast degrees of freedom and makes Eq. (15) essentially 
nonlocal in time. One may be tempted to investigate under which general conditions the 
memory kernel in (15) becomes negligible. The answer is that the full dynamics of the 
probability density p, in (9), can be properly approximated by the projected part Pp in 
(15), for those systems characterized by a rapid decay of fluctuations of the fast degrees 
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of freedom. This approximation, which is referred to in the literature as the adiabatic 
approximation [21], paves the way to determine a manifold of slow variables in the 
phase space. As already shown in previous works addressing the BE [12], the adiabatic 
approximation, which postulates a decomposition between fast and solw variables, has 
no immediate connection with the order of magnitude of the parameter e. In fact, the 
standard Chapman-Enskog technique, which relies on the adiabatic approximation in 
order to compare (5) with (6), works properly just for e 1; on the other hand, the 
method discussed in [9, 12], which also disregards the contribution of the dynamics 
of the fast variables in (15), was shown to correspond to an exact summation of the 
Chapman-Enskog expansion and does not impose any constraint on the magnitude of 
e. Now, in the spirit of the standard Chapman-Enskog technique, we may assume that 
the probability density p is a sufficiently smooth function whose time dependence is 
parameterized through a set of N phase space functions with \i = 1, N. We 

apply, then, the chain rule (summation over repeated indices is assumed) and write: 

oo r (j) 

pi>=T.^sk) 9>w (16) 

The dynamics of (A^) (provided that the boundary terms, in the integration by 
parts, are negligible) is found to be given by: 



d t (A,) = J2t k dt(A,) (17) 

where: 

d^(A^) = Jva,- rpWdT = <V^ • r)W (18) 
with T given by (10). Therefore, we obtain: 

oo 

Pp= e^dJV') (19) 

j,k=0 

with: 

a!V> = |ga!"w (20) 

It is now time to specify the various observables (A^). A crucial aspect of many 
nonequilibrium theories is, in fact, concerned with the definition of a proper set of 
coarse-grained variables triggering the dynamics of certain phenomena of interest. In 
the context of thermostatted many-particle systems, we suggest to consider the average 
phase space contraction rate (k). The latter appear to be a promising candidate, as it is 
a dynamical quantity which, in the steady state, attains a constant value related, in the 
thermodynamic limit, to the steady entropy production of Irreversible Thermodynamics 
[22]. Then, the protocol we follow is to reconstruct the dynamics of the probability 
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density by matching, at any order of e, the original microscopic dynamics given by (14) 
with the projected dynamics (19). Thus, we obtain the following equation: 



e^dj V° = - E e'r • Vp« - £ e( m+1 > (J2 • V - k) P ™ (21) 
j,fc=o 2=0 m=0 

Since p(°) evolves according to purely Hamiltonian dynamics, at the first order, e 1 , 
Eq. (21) gives: 

-r • v P « - (r ■ v - = ^(v«- r ><°> (22) 

Equation (22) addresses and defines the first order deviation from the reference 
Hamiltonian dynamics, by taking into account the effect of the external field and the 
phase space contraction induced by the thermostat. A remark is in order when discussing 
the derivative of the probability density with respect to a generic variable (A^). 
The evaluation of the average is, in principle, performed through the knowledge of the 
full probability density. This is avoided, in the Chapman-Enskog theory, by requiring 
that the collision invariants are defined only through the local Maxwellians (3). In the 
present context, instead, the assumption of weak coupling limit (e <C 1) is used to 
justify the equivalence: (A^) ~ Yli=o e *(A-*)^- ^ * s w °rth to mention that a historically 
relevant approximate solution of (22), is given by the quasi- equilibrium (or maximum- 
entropy) approximation. The technique employs a maximization of the Gibbs entropy, 
S(t) = J plogpdT, under some given constraints, and works at full equilibrium, i.e. for 
e = and in the limit t — > oo (this is, in fact, a method to derive the canonical ensemble 
in equilibrium statistical mechanics, cf. Ref. [6]). Nevertheless, it is not obvious that 
the technique also applies for e > 0, although some authors [19, 20] consider it a valid 
principle in general. Furthermore, a major drawback of this maximization method, 
outside equilibrium, is that the set of variables and the corresponding constraints 
are usually not known (typical candidates are the invariants of motion and quantities 
prescribed by the boundary conditions). We can avoid employing this approximation, 
as, due to the specific choice of variables we made, we can attempt to tackle analytically 
Eq. (22). In the limit t — > oo and for e <C 1, we seek a solution of (22) taking the 
modified Gibbs form: 

p ~ p (0) + ep (l) = e -mo(T)-e^ K (T)] (23) 

where the ip^s are Lagrange multipliers depending on the value of the chosen parameter 
(k). Let us observe that the functional form of p in (23) resembles the expression 
proposed, through a different method, also in [23]. For the projected time derivative in 
(22), one finds: 

5pW (VK ■ r )(°) ~ PK (T)(Vk ■ r ) (0) T|%y (24) 



8{k) " ' wx- u/ Q ^ 
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( A ) * ( A ) eq + e Wry (2« 



On the other hand, the computation of the microscopic time derivative, 
corresponding to the left hand side of (22), gives: 

r ■ Vp (1) = (f o ■ Vk)p4> k (25) 

Then, by inserting (24) and (25) into (22), and by integrating over the whole phase 
space, we obtain a differential equation for the ?/Vs: 

^ y <«>a) = -^[1 + 0(6)] (26) 

Therefore, up to the first order in e, the density in (23) attains the form: 

p ~ e < K > ( ' (27) 

For small external fields, the expression (27) can be further linearized in the 
nonequilibrium parameter e, which enables us to deduce the nonequilibrium steady state 
ensemble average of a generic (smooth enough) observable A for a weakly dissipative 
thermostatted particle system: 

(An) eq 

Eq. (28) is our main result, showing that the steady state average of an observable 
A comprises a first order nonequilibrium correction proportional to the equilibrium 
correlation function between phase space contraction rate and the observable itself. 
The response formula above can be interpreted as the adiabatic limit of the Green- 
Kubo formula [1], such that the fast degrees of freedom have been projected out and 
the time integral of the equilibrium time-correlation function is replaced by an effective 
equilibrium correlation function. Eq. (28) corresponds to a first order projection, on a 
manifold of slow variables, of the FDT for deterministic dissipative particle systems. 

By setting A = k and by recalling that E = <xj, with j = q(p), we used Eq. (28) 
to compute, for instance, the electrical conductivity o~(E x ) of a Gaussian thermostatted 
periodic Lorentz gas driven by an external field E x parallel to the x— axis. The numerical 
analysis performed by Lloyd et al. in Ref.[24] shows that the a(E x ) decreases nonlinearly 
with the field strength E x , which could not be justified in terms of standard linear 
response theory and was understood as a typical nonlinear (higher order) effect. In 
turn, Eq. (28) predicts, at the first order, the power law relation o ~ -7== see Fig. (1). 
Let us point out that our method applies, by construction, to dissipative systems, hence 
the external field must not vanish. Consequently, the agreement with experimental or 
numerical data is expected to be sensible in a range where the external field is finite and 
not exceedingly strong (as our result is pertinent to a first order perturbation theory). 
These expectations are, in fact, corroborated by the results shown in Fig. (2). This 
tends to suggest, a posteriori, that, in nonequilibrium response theory for thermostatted 
particle systems, the adiabatic limit stands as a plausible physical approximation. 
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Figure 1. Electrical conductivity a vs. electric field E x . 

Black dots: numerical results of Lloyd et al. [24]. 

Blue line: adiabatic approximation of the Liouville Equation 
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0.85 0.90 



Figure 2. Magnification of Fig. (1) in the range where the theory 
best fits the numerical data. 



4. Conclusions 

We remark that the novelty of the proposed approach stems from a suitable 
combination of a perturbation theory with a projection operator technique in the 
context of deterministically thermostatted particle systems with nonvanishing phase 
space contraction rate. Our work aimed to show that the hypothesis of a time scales 
separation in a thermostatted particle system enables to express the expectation value 
of an observable solely in terms of few selected relevant variables (in our case the average 
phase space contraction rate) which drive the dynamics of the probability density. The 
issue of determining response formulae in terms of few relevant observables is crucial 
to establish a bridge between a macroscopic hydro dynamic-like description and the 
underlying Liouville Equation. Alternative derivations which lead to response formulae 
deduced from the exact Liouville propagator, without resorting to any projection 
technique, are known [25]. We believe it is convenient to also look at approximate 
formulae which might offer a complementary point of view and might be, in some cases, 
of more immediate application, as we showed above in computing the conductivity 
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of the Lorentz gas. It is also worth to notice that, regardless of the intensity of the 
external field, for any finite field the steady state attractor has a lower dimension than 
the embedding space. Thus, since the support of the steady state measure lies on the 
attractor, the measure is not smooth and one has to resort to SRB measures [4]. On 
the other hand, the procedure of projecting out the fast degrees of freedom, is believed 
to ensure regularity to the resulting density, as also discussed in [26]. Therefore, the 
suggested technique seems, also, to corroborate the applicability of the method proposed 
by Vulpiani et al. in [2], without invoking unstable manifolds of Anosov systems or 
introducing any noise occurring in not isolated physical systems 
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